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Abstract.  Bentley,  Weide  and  Yao  have  proposed  an  expected  0(N ) cell  tech- 
nique for  computing  Voronoi  diagrams  in  two  dimensions  that  does  not  generalize 
readily  to  three.  In  this  paper  their  work  is  further  developed  and  generalized 
to  produce  expected  0(N ) and  0(N 4/3)  algorithms  for  constructing  Voronoi  di- 
agrams in  two  and  three  dimensions,  respectively.  Computational  experience  is 
presented  for  the  algorithm  in  two  dimensions. 
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1.  Introduction 

Consider  a finite  set  5 of  points  in  Euclidean  space,  and  let  P be  any  point 
in  5.  The  Voronoi  polyhedron  V(P,S)  of  P relative  to  S is  the  set  of  all 
points  in  the  space  such  that  P is  as  close  to  any  point  in  this  set  as  is 
any  other  point  in  S.  V(P,  S)  is  the  result  of  a bisection  process  for  P 
with  respect  to  S;  i.  e.,  V[P,S)  is  the  intersection  of  the  closed  half-spaces 
that  contain  P and  that  are  determined  by  the  perpendicular  bisectors  of 
the  line  segments  connecting  P with  the  other  points  in  5.  It  follows  that 
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V(P,  S)  is  a closed  convex  polyhedron,  possibly  unbounded,  which  contains 
P in  its  interior.  The  Voronoi  polyhedra  of  the  points  in  S (relative  to  5 if 
not  stated  otherwise)  fill  the  space  without  common  interior  points,  and  the 
union  of  their  boundaries  forms  a diagram,  the  Voronoi  diagram  for  5,  which 
partitions  the  space  into  the  Voronoi  polyhedra.  A point  in  the  space  that 
is  a vertex  of  some  Voronoi  polyhedron  is  said  to  be  a vertex  of  the  Voronoi 
diagram  for  5.  In  ^-dimensional  Euclidean  space  a vertex  of  the  Voronoi 
diagram  for  a set  S is  always  the  common  vertex  of  at  least  k + 1 Voronoi 
polyhedra  of  points  in  5.  It  is  called  a degenerate  vertex  whenever  it  is  a 
vertex  of  more  than  k + 1 Voronoi  polyhedra. 

Bentley,  Weide  and  Yao  [l]  have  shown  how  the  Voronoi  diagram  for  a 
set  5 of  N points,  chosen  independently  from  a uniform  distribution  on  a 
square  in  the  plane,  can  be  constructed  in  linear  expected  time.  With  M 
defined  as  the  largest  integer  less  than  or  equal  to  V1/2,  i.  e.  the  floor  of  TV1/2, 
Bentley,  et  al.  divide  the  square  into  M2  equal-sized  square  cells.  With  L(N ) 
defined  as  the  floor  of  log  N (here  log  denotes  the  natural  logarithm),  cells  not 
contained  in  the  outermost  L(N ) layers  of  cells  of  the  square  are  called  ‘inner 
cells’,  the  rest  ‘outer  cells’.  Essentially,  the  algorithm  takes  three  steps.  In 
the  first  step,  each  point  in  S is  assigned  (in  constant  time)  to  a cell  in  which 
it  is  contained.  In  the  second  step,  the  Voronoi  polygons  of  points  assigned 
to  inner  cells  are  constructed.  Given  a point  P assigned  to  an  inner  cell, 
a search  is  conducted  through  each  of  the  layers  of  cells  surrounding  P for 
points  assigned  to  cells  in  these  layers.  This  search  procedure,  called  a spiral 
search  around  P,  starts  with  the  cell  assigned  to  P,  and  then  proceeds  in 
outward  direction  to  each  of  the  layers  of  cells  surrounding  this  cell.  Vr(P,  5) 
is  progressively  built  through  a bisection  process  for  P with  respect  to  the 
set  of  points  in  5 found  through  the  search.  A geometric  test  is  available 
for  deciding  whether  V(P,  5)  has  been  obtained.  In  most  cases  Vr(P,  5)  is 
obtained  after  examining  only  a small  number  of  cells  and  points  so  that  the 
expected  time  for  the  second  step  is  O(N).  Finally,  in  the  third  step,  the 
Voronoi  polygons  of  those  points  assigned  to  outer  cells  are  constructed  by 
applying  to  them  an  0(K  log  K)  worst  case  algorithm,  e.  g.  Shamos’  [3]. 
The  expected  time  for  this  last  step  is  also  O(JV). 

In  this  paper,  we  present  a version  of  this  algorithm  in  which  the  third 
step  is  replaced  by  steps  that  also  require  expected  time  O(N)  but  that  have 
the  advantage  of  generalizing  to  three  dimensions.  Again,  with  M defined  as 
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the  floor  of  TV1/2,  the  square  is  divided  into  M2  equal-sized  cells.  We  define 
inner  cells  and  outer  cells  differently  from  [ 1] . Namely,  cells  contained  in  the 
outermost  two  layers  of  cells  of  the  square  are  now  called  outer  cells , the 
rest  inner  cells.  As  in  [l],  each  point  in  5 is  assigned  to  a cell  in  which  it  is 
contained.  Voronoi  polygons  of  points  assigned  to  inner  cells  are  constructed 
by  applying  the  Bentley,  et  al.  cell-based  method  together  with  a generaliza- 
tion of  its  geometrical  test.  Finally,  a modified  version  of  Bowyer’s  insertion 
process  [2]  is  used  to  handle  points  assigned  to  outer  cells. 

An  insertion  process  is  a method  for  updating  the  Voronoi  diagram  for  a 
set  of  k — 1 points  for  an  additional  kth  point  in  order  to  obtain  the  Voronoi 
diagram  for  all  k points.  The  first  step  in  Bowyers  algorithm  consists  of 
identifying  a vertex  of  the  Voronoi  diagram  for  the  k — 1 points  that  will  not 
be  a vertex  of  the  Voronoi  diagram  for  the  k points.  To  do  this,  Bowyer  uses 
a ‘walk’  that  starts  near  the  centroid  of  the  k — 1 points  and  that  ends  at  a 
point  in  whose  Voronoi  polygon  the  kth  point  lies.  The  point  thus  found  is 
clearly  a point  among  the  k — 1 points  that  is  closest  to  the  kth  point.  We 
have  modified  this  step  by  taking  advantage  of  the  cell  structure.  In  place  of 
a walk  we  use  a spiral  search  around  the  kth  point  and,  as  described  in  [l], 
find  a point  among  the  k — 1 points  that  is  closest  to  it.  The  remainder  of 
Bowyer’s  algorithm,  in  which  vertices  of  the  Voronoi  diagram  for  the  k — 1 
points  which  lie  in  the  Voronoi  polygon  of  the  kth  point  are  deleted,  and  the 
vertices  of  the  Voronoi  polygon  of  the  kth  point  are  added,  is  left  essentially 
unaltered. 

In  Section  2,  the  two-dimensional  algorithm  is  outlined  and  justified.  In 
Section  3,  we  prove  the  algorithm  is  of  expected  O(N)  complexity.  In  Sec- 
tion 4,  the  algorithm  is  generalized  to  three  dimensions  and  shown  to  be  of 
expected  0(TV4//3)  complexity.  Finally,  in  Section  5,  computational  experi- 
ence with  the  implementation  of  the  two-dimensional  algorithm  is  presented. 


2.  The  two-dimensional  algorithm 

In  the  following,  cells  are  divided  into  four  classes  of  cells  (Figure  l).  Inner 
cells  not  contained  in  any  of  the  outermost  L(N ) layers  of  cells  of  the  square 
are  called  class  1 cells.  Inner  cells  within  T(TV)  layers  of  cells  from  exactly 
one  side  of  the  square  are  called  class  2 cells.  Inner  cells  within  T(TV)  layers 
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Figure  1:  Regions  of  square  that  contain  the  four  classes  of  cells.  Class  1,  2, 
and  3 cells  are  inner  cells,  and  class  4 cells  are  outer  cells. 
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of  cells  from  two  sides  of  the  square  are  called  class  3 cells.  Finally,  all  outer 
cells  are  called  class  4 cells.  We  assume  that  N is  large  enough  so  that  none 
of  the  four  classes  is  empty.  We  define  54  C 5 as  the  set  of  points  assigned 
to  class  1 cells.  52,  53,  54  are  analogously  defined  with  respect  to  class  2, 
class  3,  class  4 cells,  respectively. 

Throughout  the  outline  of  the  algorithm  that  follows,  we  say  a point  in 
5 has  been  processed  if  its  Voronoi  polygon  has  been  constructed.  We  say 
a cell  is  empty  if  no  points  of  5 are  assigned  to  it.  We  say  an  inner  cell 
has  been  activated  if  it  has  been  found  to  be  empty  or  if  each  point  of  5 
assigned  to  the  cell  has  been  processed.  We  select  a centermost  class  1 cell. 
The  order  in  which  inner  cells  are  activated  by  our  algorithm  is  determined 
by  proceeding  through  the  layers  of  cells  in  the  square  in  a spiral-like  fashion 
around  the  centermost  cell.  We  say  an  inner  cell  is  the  currently  active  cell  if 
the  points  assigned  to  it  are  currently  being  processed.  We  say  a point  in  the 
currently  active  cell  is  the  currently  active  point  if  the  point  is  currently  being 
processed.  We  say  a cell  has  been  searched  during  a spiral  search  around  a 
point  in  5 if  all  points  assigned  to  the  cell  have  already  been  found  through 
the  search. 

We  say  that  points  P and  P'  in  5 are  Voronoi  neighbors  in  S if  V(P,  5) 
and  V(P',S)  are  contiguous  in  the  Voronoi  diagram  for  5.  For  each  P in 
5 we  let  Vr(P,  5)  represent  the  set  of  Voronoi  neighbors  of  P in  5,  and  we 
assume  V(P,  5)  can  be  readily  extracted  from  V(P,  5).  During  the  execution 
of  the  algorithm,  for  each  P in  5 that  has  not  been  processed  we  define  ZV(P) 
as  the  set  of  points  assigned  to  inner  cells  that  are  known  Voronoi  neighbors 
of  P in  5;  i.  e.,  the  set  of  points  in  5\S4  that  have  been  processed  and  found 
to  be  Voronoi  neighbors  of  P in  5.  Every  point  assigned  to  an  inner  cell  is 
processed  by  our  algorithm  with  a spiral  search  around  the  point.  Given 
P,  one  such  point,  and  assuming  it  is  the  currently  active  point,  we  let  dt 
represent  the  radius  of  the  largest  circle  centered  at  P whose  interior  only 
intersects  cells  that  have  been  searched  up  to  the  current  time.  We  also  let  Z* 
be  the  set  of  points  in  5 that  have  been  found  up  to  the  current  time  during 
the  search.  Then,  as  in  [1],  we  define  the  tentative  Voronoi  polygon  V1  of  P 
at  the  current  time  as  the  Voronoi  polygon  of  P relative  to  {P}  U Zl  U Z*(P) 
(the  plane  if  Zl  and  ZV(P)  are  empty).  Accordingly,  when  the  search  starts 
and  Zl  is  still  empty,  we  define  the  initial  tentative  Voronoi  polygon  V1  of  P 
as  the  tentative  Voronoi  polygon  of  P at  that  time.  Also  as  in  [l],  each  Vl  is 
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constructed  by  updating  a previous  one  in  such  a way  that  each  V1  is  exactly 
the  result  of  a bisection  process  for  P with  respect  to  Zl  U Zv  (P)- 

Given  points  0 and  O'  in  the  plane,  we  let  QQ ' and  dist(Q ,0')  repre- 
sent, respectively,  the  closed  linear  segment  that  connects  0 and  Q',  and  the 
distance  from  0 to  O'.  In  the  outline  of  the  algorithm,  given  Q,  a point  in 
the  plane,  and  W , a finite  or  infinite  subset  of  the  plane,  we  let  dmax(Q , W) 
represent  the  maximum  value  of  dist(Q ,Q')  for  Q'  in  W.  Given  rays  rfi  and 
f2  of  common  origin,  we  let  771(7*1,7*2)  represent  the  size  of  the  angle  produced 
by  a clockwise  rotation  from  rfi  to  f2.  Given  points  Q and  O’  in  the  plane,  we 
let  QQ ' represent  the  ray  through  O'  of  origin  0.  Given  a point  P in  5 2,  we 
let  both  1\{P)  and  12{P ) represent  the  ray  of  origin  P that  is  perpendicular 
to  the  side  of  the  square  closest  to  P.  Given  a point  P in  S3,  we  let  li(P) 
and  l2 (P)  represent  the  two  rays  of  origin  P that  are  perpendicular  to  the 
two  sides  of  the  square  closest  to  P with  m(/1(P),  l2(P))  equal  to  90°.  Given 
P,  a point  in  S2  or  S3,  and  assuming  it  is  the  currently  active  point,  we 
let  mti  and  mt2  represent  at  the  current  time  the  smallest  positive  values 
of  m(PP',  l\(P))  and  m(/2(P),  PP'),  respectively,  for  P'  in  Zl  U ZV{P).  We 
say  points  Pi  and  P2  in  Zl  U ZV(P)  determine  mtl  and  mt2 , respectively,  if 
mn  equals  m(P P1.l1(P))  and  m<2  equals  m(/2(P),  P P2).  Assuming  points 
Pi  and  P2  determine  rritl  and  m(2,  respectively,  we  let  Cl  represent  the  inte- 
rior of  that  region  of  the  plane  obtained  by  a clockwise  rotation  from  P Pi  to 
PP2.  We  are  now  ready  to  formulate  the  algorithm. 

Start  of  algorithm. 

Step  1.  Assign  points  to  cells  and  select  first  class  1 cell  to  be  activated. 

Let  M be  the  floor  of  N1^2. 

Partition  the  square  into  M 2 equal-sized  square  cells. 

Determine  inner  and  outer  cells. 

Assign  each  point  in  5 to  a cell. 

For  each  cell,  list  the  points  assigned  to  it. 

Determine  the  centermost  cell. 

If  the  centermost  cell  is  empty  then  go  to  Step  2. 

Else  designate  this  cell  as  the  currently  active  cell. 

Go  to  Step  3. 

Step  2.  Select  next  inner  (class  1,  2,  or  3)  cell. 
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If  all  inner  cells  have  been  activated  then  go  to  Step  8. 

Else  choose  the  next  inner  cell  to  be  activated. 

If  this  cell  is  empty  then  go  to  Step  2. 

Else  designate  this  cell  as  the  currently  active  cell. 

Determine  class  for  currently  active  cell. 

If  class  1 then  go  to  Step  3. 

Else  go  to  Step  4. 

Step  3.  Construct  Voronoi  polygon  of  a point  in  S i. 

Let  P be  a point  assigned  to  the  currently  active  cell  that  has  not  been 
processed. 

Designate  P as  the  currently  active  point. 

Start  spiral  search  around  P and  construct  V1. 

Update  V1  and  each  subsequent  V1  as  appropriate. 

For  each  \n  compute  Dt  = dmax(P,Vt ) and  dt. 

Terminate  search  when  one  of  the  following  criteria  is  met. 

1.  2 Dt  < dt . 

2.  All  cells  in  the  square  have  been  searched. 

Upon  termination  go  to  Step  7. 

Step  4.  Begin  construction  of  Voronoi  polygon  of  a point  in  S 2 or  S 3. 

Let  P be  a point  assigned  to  the  currently  active  cell  that  has  not  been 
processed. 

Designate  P as  the  currently  active  point. 

Determine  l\(P)  and  Z2(P). 

Start  spiral  search  around  P and  construct  V1. 

Update  Ul  and  each  subsequent  Vf  as  appropriate. 

For  each  Uf  compute  Dt  = dmax{P.Vt).  dtl  mfl,  and  77if2,  and  deter- 
mine CL 

Terminate  search  when  one  of  the  following  criteria  is  met. 

1.  2 Dt  < dt . 

2.  All  cells  in  the  square  have  been  searched. 

3.  All  cells  that  intersect  Cl  have  been  searched. 


7 


Upon  termination,  if  neither  criterion  1 nor  criterion  2 has  been  met 
then  go  to  Step  5. 

Else  go  to  Step  7. 

Step  5.  Continue  construction  of  Voronoi  polygon  of  a point  in  S 2 or  S3. 
Let  C = CV 

Let  Pi  and  P2  be  points  in  S that  determine  mtl  and  mt 2,  respectively. 
Compute  d!  = dmax(P,  {Pi,  i^})- 
Resume  spiral  search  around  P. 

Update  each  Vrt  as  appropriate. 

For  each  V1  compute  Dt  = dmax(P,Vt ) and  dt. 

Terminate  search  when  one  of  the  following  criteria  is  met. 

1.  2 Dt  < dt . 

2.  All  cells  in  the  square  have  been  searched. 

3.  d!  < dt 

Upon  termination,  if  neither  criterion  1 nor  criterion  2 has  been  met 
then  go  to  Step  6. 

Else  go  to  Step  7. 

Step  6.  Complete  construction  of  Voronoi  polygon  of  a point  in  S 2 or  S3. 
Resume  spiral  search  around  P . 

Update  each  V1  as  appropriate. 

For  each  \n  compute  Dt  = dmax(P , V1  \ C)  and  dt. 

Terminate  search  when  one  of  the  following  criteria  is  met. 

1.  2 Dt  < dt. 

2.  All  cells  in  the  square  have  been  searched. 

Upon  termination  go  to  Step  7. 

Step  7.  Save  Voronoi  polygon  of  a point  assigned  to  an  inner  cell. 

Identify  U(P,  S)  with  V1 . 

Mark  P as  processed  and  save  V(P,  S). 

For  each  P'  in  V(P,  S)  that  has  not  been  processed  let 
ZV(P ')  = ZV(P')  U {P}. 
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Determine  whether  currently  active  cell  has  been  activated. 

If  activated  then  go  to  Step  2. 

Else  if  P is  in  Si  then  go  to  Step  3. 

Else  go  to  Step  4. 

Step  8.  Construct  and  save  Voronoi  polygons  of  points  in  S 4. 

Determine  S4. 

If  S4  is  empty  then  stop. 

Else  perform  insertion  process  on  S4. 

Perform  for  each  P in  S4  a bisection  process  with  respect  to 

V(P,  S4)  U ZV(P)  and  identify  V(P,  S)  with  result  of  process. 

For  each  P in  S4  mark  P as  processed  and  save  V(P,  5). 

Stop. 

End  of  algorithm. 

Justification  of  algorithm.  As  established  in  [l],  the  Voronoi  polygons  of 
points  in  Si  can  be  constructed  with  Step  3 of  the  algorithm.  Let  P,  V4,  Z)t, 
dt  be  as  defined  in  Step  3.  Let  P'  be  a point  in  5 with  2 Dt  < dist(P , P').  It 
follows  that  P'  can  not  affect  V1  since  the  perpendicular  bisector  of  P P'  does 
not  intersect  V1.  Thus,  during  the  spiral  search  around  P,  we  may  conclude 
that  V(P,  S ) is  equal  to  V1  as  soon  as  2 Dt  < dt. 

We  justify,  with  the  aid  of  Figure  2,  that  we  can  produce  the  Voronoi 
polygons  of  points  in  S2  or  S3  with  Steps  4,  5,  6 of  the  algorithm.  Let  P be 
as  defined  in  Step  4,  and  let  C,  Pi,  P2,  d!  be  as  defined  in  Step  5.  Let  V 
be  the  portion  of  V(P,  {P,  Pi,  P2})  that  is  contained  in  C (shaded  region  in 
Figure  2).  Let  P'  be  a point  in  S that  does  not  lie  in  C or  in  the  interior  of 
the  circles  with  centers  (P  + Pi)/2,  (P  + P2)/2,  and  diameters  dist(P1  Px), 
dist(P,  P2),  respectively,  as  shown  in  Figure  2.  It  follows  that  P'  does  not 
affect  V as  indicated  in  Figure  2 by  the  perpendicular  bisector  b of  PPL 
We  note  that  points  in  S contained  in  C are  found  through  the  spiral  search 
around  P with  Step  4,  and  that  points  in  S contained  in  the  circle  with  center 
P and  radius  d!  are  found  with  Step  5 (this  circle  contains  the  two  circles 
mentioned  above  and  is  easier  to  search).  Thus,  with  Step  6,  in  which  V is 
not  considered  during  the  geometrical  test  and  through  which  only  points 
that  do  not  affect  V are  found,  we  can  produce  V(P,  5). 
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Figure  2:  P\  a point  in  5 that  does  not  affect  V (shaded  region). 
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Finally,  we  can  construct  the  Voronoi  polygons  of  points  in  S4  with  Step  8 
of  the  algorithm  since  two  points  in  S4  that  are  Voronoi  neighbors  in  5 must 
be  Voronoi  neighbors  in  S4. 


3.  Proof  of  complexity 

In  this  section,  we  assume  that  the  points  in  5 have  been  chosen  from  a 
uniform  distribution  on  the  square,  and  prove  that  the  algorithm  presented 
in  Section  2 has  expected  O(N)  execution  time.  First,  we  state  the  following 
theorem  by  Bentley,  Weide  and  Yao  [1].  Here,  the  time  involved  is  defined 
as  the  number  of  cells  plus  the  number  of  points  examined. 

Theorem  Let  P be  a point  in  the  square,  and  let  P'  be  a point  in  S closest  to 
P . Then  the  expected  time  required  to  find  P1  through  a spiral  search  around 
P is  constant. 

By  modifying  the  proof  of  the  theorem,  Bentley,  et  al.  also  establish  the 
following  two  observations  crucial  to  our  proof  of  optimality. 

1 . Let  P be  a point  in  S such  that  at  least  one  point  in  S is  contained  in  each 
of  the  octants  around  P shown  in  Figure  3.  Then  the  expected  time  required 
to  find  at  least  one  point  in  each  octant  through  a spiral  search  around  P is 
constant. 

2.  Under  similar  assumptions,  let  Pl}  i = 1 , . . . , 8,  be  the  first  eight  points 

in  S obtained  through  a spiral  search  around  P such  that  they  are  con- 
tained in  the  octants  around  P , one  per  octant.  Let  d and  D be  the  values 
of  dmax(P,{Pi, ...,  P»})  and  dmax(P.  V(P,  {P,  . . . , P8}));  respectively. 

V(P,  S)  can  be  constructed  by  searching  only  those  cells  intersecting  the  inte- 
rior of  the  circle  with  center  P and  radius  2 D (Figure  3).  Since  2 D <y/2d, 
it  follows  from  Observation  1 that  the  expected  time  required  to  search  all  of 
these  cells  through  a spiral  search  around  P is  constant. 

Assigning  the  points  in  5 to  the  appropriate  cells  can  be  accomplished  in 
0(N ) time  [l].  Thus,  it  will  suffice  to  show  that  the  expected  time  involved  in 
constructing  wdth  the  algorithm  all  Voronoi  polygons  of  points  in  S;,  for  each 
i,  i = 1, . . . , 4,  is  bounded  above  by  O(N).  Let  P be  a point  in  5 \ S4  so  that 
V(P,  S)  is  constructed  with  the  algorithm  through  a spiral  search  around  P. 
In  what  follows,  we  let  w denote  the  time  involved  in  constructing  V(P,  5), 
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Figure  3:  Points  Pt,  i = contained  in  octants  around  P.  P'  is 

outside  the  circle  of  radius  2 D so  it  does  not  affect  V[P,  {P,  Px, . . . , P&}) 
as  indicated  by  the  perpendicular  bisector  b of  PP'. 
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and  use  the  fact  that  w is  bounded  above  by 

o(j  + 

i=i 

where  j is  the  number  of  cells  examined  with  the  search,  k is  the  number 
of  points  in  those  j cells,  and  Vi  is  the  number  of  vertices  in  the  tentative 
Voronoi  polygon  of  P when  the  ith  point  is  found  through  the  search.  Finally, 
we  let  E(w)  denote  the  expected  value  of  w,  i.  e.  the  expected  time  involved 
in  constructing  Vr(P,  S). 

Proof  for  S\.  Let  P be  a point  in  Si.  As  in  [1],  we  say  P is  closed  if  at 
least  one  point  in  5 is  contained  in  each  of  the  octants  around  P as  shown 
in  Figure  3. 

Let  pi  be  the  probability  that  P is  closed,  and  t\  the  expected  number  of 
points  examined  while  constructing  V(P,  S)  with  the  algorithm  when  P is 
closed.  P2  and  t2  are  similarly  defined,  respectively,  for  P not  closed. 

If  t is  the  expected  number  of  points  examined  while  constructing  V(P,  5) 
with  the  algorithm,  then  t = pi  • t\  + p2  * ^2- 

Of  course  pi  < 1,  and  from  Observation  2 above,  ti  = 0(1).  Since  at  most 
all  points  are  examined  when  P is  not  closed,  it  follows  that  t2  < 0(N). 
Next,  we  find  an  upper  bound  for  p2  using  an  argument  of  [1].  If  no  points 
are  found  in  a given  octant,  then  at  least  0(L(N)2)  cells  must  be  empty.  The 
probability  of  the  octant  being  empty  is  then  bounded  above  by  \ 

It  follows  that  p2  < 8 
Therefore, 

t = 1 • 0(1)  + 8e~0(t(Ar)2)  • O(N)  = 0(1). 

A similar  argument  can  be  used  to  show  that  the  expected  number  of  cells 
examined  while  constructing  V(P,  5)  with  the  algorithm  is  constant.  Finally, 
since  the  number  of  vertices  in  any  tentative  Voronoi  polygon  of  P is  at  most 
the  number  of  points  examined  while  constructing  V(P,  5),  the  expected 
number  of  vertices  in  any  tentative  Voronoi  polygon  of  P is  also  constant. 

It  follows,  then,  that 


E{w)  < 0(1)  + 0(1)  • 0(1)  = 0(1) 

for  each  point  in  Si. 

Since  at  most  N points  are  contained  in  5 1 then  the  expected  time  required 
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for  Si  is 


N • 0(1)  = O(N). 

Proof  for  S2.  Let  P be  a point  in  S2  and,  without  any  loss  of  generality, 
assume  P is  within  L(N)  layers  of  cells  from  the  right-hand  side  of  the 
square.  Let  l represent  the  ray  that  both  h(P)  and  Z2(P)  represent.  As 
shown  in  Figure  4,  we  say  P is  closed  if  within  the  first  L{N ) layers  of  cells 
that  surround  P at  least  one  point  in  S is  contained  in  each  of  six  octants 
around  P,  octants  I through  V/,  and  at  least  one  point  in  5 (which  may  be 
one  of  the  points  in  octants  I or  VI)  is  found  in  each  of  the  upper  and  lower 
portions  of  the  outermost  layer  of  the  square. 

If  P is  closed  let  Pi  and  P2  be  points  in  5 within  the  first  L(N ) layers  of  cells 

that  surround  P in  the  upper  and  lower  portions  of  the  outermost  layer  of 

— * — ♦ 

the  square,  respectively,  with  the  smallest  positive  values  of  m(PPi,l)  and 
m(/,  PP2).  Let  C be  the  interior  of  the  region  of  the  plane  obtained  by  a 
clockwise  rotation  from  PPi  to  PP2. 

Clearly,  since  P is  further  than  two  cells  from  all  sides  of  the  square,  PPi  and 
PP2  intersect  the  boundary  of  the  square  at  points  within  t'he  first  2 L(N) 
layers  of  cells  that  surround  P.  Therefore,  by  examining  the  first  2 L(N) 
layers  of  cells  that  surround  P,  all  cells  intersecting  C are  examined. 

Let  d'  be  the  value  of  dmax(P,  {Pi,  P2}).  Then  the  circle  of  radius  d'  with 
center  P is  also  contained  in  the  first  2 L(N)  layers  of  cells  that  surround  P. 
Finally,  let  U be  the  Voronoi  polygon  of  P relative  to  those  points  within 
the  first  L(N)  layers  of  cells  that  surround  P.  Note  that  points  outside  C 
and  the  circle  of  radius  d'  with  center  P do  not  affect  the  part  of  U that  is 
contained  in  C.  Accordingly,  let  D'  be  the  value  of  dmax(P,U  \ C).  Then, 
since  P is  closed,  the  circle  of  radius  2 D'  with  center  P is  also  contained  in 
the  first  2 L(N)  layers  of  cells  that  surround  P. 

It  follows  from  these  observations  that  the  Voronoi  polygon  of  a closed  point 
can  be  constructed  with  the  algorithm  by  examining  no  more  than  the  first 
2 L(N)  layers  of  cells  that  surround  the  point. 

Let  pi  be  the  probability  that  P is  closed,  and  ti  the  expected  number  of 
points  examined  while  constructing  V(P,  5)  with  the  algorithm  when  P is 
closed.  p2  and  t2  are  similarly  defined,  respectively,  for  P not  closed. 

If  t is  the  expected  number  of  points  examined  while  constructing  V(P,S ) 
with  the  algorithm,  then  t = pi  • ti  -f  p2  • t2. 


14 


Outermost 
2 layers 


L(N)  cells 


/" 


pr 


\ 


\ 


Octant  II 


J2 

o 

o 


\ 


\ 


\ # 

\ 

\ 


\ 


Octant  III 


\ 


\ 


\ 


\ 


/ 


Octant  IV  / 

/ 

/ 

/ 

/ 

/ 

/ Octant  V 

/ 


/ 


Id 


eta 


nt  I 


/ 


/ 


/ 


/ 


/ 


/ 


Clipper  portion 


i 


\ 


\ 


\ 


\ 


\ 


> Lower  portion 

\ 


\ 


eta 

VI 


nt 


\ 


Figure  4:  A closed  point  in  S2.  One  point  is  contained  in  each  of  the  octants 
I through  VI.  Pi  and  P2  are  contained  in  the  upper  and  lower  portions  of 
the  outermost  layer  of  the  square,  respectively. 
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Of  course  pi  < 1,  t2  < O(N),  and  from  the  above  discussion  0(L(N)2)  is  an 
upper  bound  for  t j. 

In  order  to  find  an  upper  bound  for  p2,  we  argue  as  follows.  If  no  points  are 
found  in  one  of  the  octants  II  through  V,  then  at  least  0(L(N)2)  cells  must 
be  empty.  If  no  points  are  found  in  one  of  the  octants  I and  VI,  then  at 
least  0(L(N ))  cells  must  be  empty.  Finally,  if  either  the  upper  or  the  lower 
portion  of  the  outermost  layer  of  the  square  is  empty,  then  0(L(N ))  cells 
must  be  empty.  It  follows  that  p2  1 4 e-o(l'(-^)2)  4.  4 e-°{LiN))m 
Thus, 

t < 1 • 0{L(Nf)  + (4e-0{L<-N),)  + 4e~°WN)))  ■ O(N) 

= 0(L(N)2)  + 0(1)  = 0{L(Nf). 

An  argument  similar  to  that  used  for  points  in  S\  can  be  used  to  show  that 
the  expected  number  of  cells  examined  while  constructing  V(P,  S)  with  the 
algorithm  is  0{L{N)2),  and  that  the  expected  number  of  vertices  in  any 
tentative  Voronoi  polygon  of  P is  also  0(L(N)2).  Thus, 

E(w)  < 0(L(N)2)  + 0(L(N)2)  ■ 0(L(N)2)  = 0(L(N)4) 

for  each  point  in  S2. 

The  expected  number  of  points  in  S2  is  0{Nl^2L[ N)).  Hence,  the  expected 
time  required  for  S2  is 

0(N1/2l(N))  ■ 0(L(N )4)  = O(N). 

Proof  for  S3.  Let  P be  a point  in  S3  and,  without  any  loss  of  generality, 
assume  P is  within  L(N)  layers  of  cells  from  the  right-hand  and  bottom 
sides  of  the  square.  Let  /j  and  /2  represent  the  rays  that  l\(P)  and  /2(P) 
represent,  respectively.  As  shown  in  Figure  5,  we  say  P is  closed  if  within 
the  first  L(N)  layers  of  cells  that  surround  P at  least  one  point  in  S is 
contained  in  each  of  four  octants  around  P,  octants  / through  IV,  and  the 
right-hand  and  bottom  portions  of  the  outermost  layer  of  the  square. 

That  O(N)  is  an  upper  bound  for  the  expected  time  required  for  S3  now 
follows  by  an  argument  similar  to  the  one  used  for  S2. 

Proof  for  S 4.  The  expected  time  required  by  the  insertion  process  of  Step  8 of 
the  algorithm  is  at  most  proportional  to  the  product  of  the  expected  number 
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Figure  5:  A closed  point  in  5 3.  One  point  is  contained  in  each  of  the  octants  / 
through  I\  . Pi  and  Pn  are  contained  in  the  right-hand  and  bottom  portions 
of  the  outermost  layer  of  the  square,  respectively. 
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of  points  in  54  and  the  expected  maximum  number  of  vertices  in  the  Voronoi 
diagram  for  any  subset  of  S4.  Since  the  expected  number  of  points  in  S4  is 
OfiV1/2),  it  follows  from  the  Euler-Poincare  formula  that  this  time  is  at  most 


0(N1/2)  ■ 0(N1/2)  = O(N). 


Finally,  it  suffices  to  show  that  the  expected  time  required  by  the  bisection 
process  of  Step  8 is  also  at  most  O(N).  Let  r be  the  number  of  points  in  S4. 
Let  i = 1, ...  ,r,  be  the  points  in  S4.  For  each  i,  i = 1 , . . . ,r,  define  as 
the  number  of  points  in  V’(Pn  S4),  Ui  as  the  final  number  of  points  in  ZV(P;), 
and  Vi  as  the  maximum  number  of  vertices  of  any  polygon  obtained  during 
the  bisection  process  for  Pt. 

It  follows  that  the  time  required  by  the  bisection  process  is  at  most  propor- 
tional to 

r r r 

£K-  + Ui)  • Vi  < £(u>*  + Ui)2  = £(w?  + 2wiUi  + u2). 

i = l i= 1 i=l 

Again,  r has  expected  value  0(N1'2),  so  that  by  the  Euler-Poincare  formula, 
J2i=iwi  ^as  expected  value  0(iV1//2). 

In  order  to  calculate  an  upper  bound  for  the  expected  value  of  each  u j, 
i = 1, . . . ,r,  we  proceed  as  follows.  Given  z,  1 < i<  r,  let  P'  be  a point  in 
5 not  contained  in  54  such  that  P'  is  outside  the  first  2 L(N)  layers  of  cells 
that  surround  P;.  As  previously  proven,  P'  is  a Voronoi  neighbor  of  Pl  in 
S with  probability  at  most  proportional  to  so  that  the  expected 

value  for  Ui  is  bounded  above  by 

0(L(N)2)  • 1 + N-  e-0[L{N))  = 0(L(N)2)  + 0(1)  = 0(L(N)2). 

It  follows  now  that  the  expected  value  of  Ya=i{w{  + 2 WiUi  -f  u\)  is  at  most 

0{N)  + 0{N1/2)  • 0(L(N)2)  + 0{ N1/2)  • 0{L{N )4)  = 0(N); 

•• 

i.  e.,  the  expected  time  required  by  the  bisection  process  is  at  most  0(N). 
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4.  The  three-dimensional  algorithm 

We  now  present  an  algorithm  for  constructing  the  Voronoi  diagram  for  a 
set  5 of  N points  contained  in  a cube  in  three-dimensional  Euclidean  space. 
First,  with  M defined  as  the  floor  of  N 1//3  we  divide  the  cube  into  M 3 equal- 
sized cubic  cells.  Cells  contained  in  the  outermost  two  layers  of  cells  of  the 
cube  we  call  outer  cells , the  rest  inner  cells.  Next,  as  in  the  two-dimensional 
case,  each  point  in  5 is  assigned  to  a cell  in  which  it  is  contained.  Finally, 
the  Voronoi  polyhedron  of  each  point  in  S is  constructed  according  to  its  cell 
assignment  by  generalizing  the  two-dimensional  algorithm  of  Section  2. 

In  order  to  outline  the  algorithm,  cells  are  further  divided  into  five  classes 
of  cells.  With  L(N ) defined  again  as  the  floor  of  log  TV,  inner  cells  not  con- 
tained in  any  of  the  outermost  L(N)  layers  of  cells  of  the  cube  are  called 
class  1 cells.  Inner  cells  within  L(N ) layers  of  cells  from  exactly  one  face  of 
the  cube  are  called  class  2 cells.  Inner  cells  within  L(N)  layers  of  cells  from 
exactly  two  faces  of  the  cube  ared  called  class  3 cells.  Inner  cells  within  L(N ) 
layers  of  cells  from  three  faces  of  the  cube  are  called  class  4 cells.  Finally,  all 
outer  cells  are  called  class  5 cells.  We  assume  N is  large  enough  so  that  none 
of  the  five  classes  is  empty.  We  define  Si  C 5 as  the  set  of  points  assigned  to 
class  1 cells.  5 2,  S3,  S4,  S5  are  analogously  defined  with  respect  to  class  2, 
class  3,  class  4,  class  5,  respectively. 

Throughout  the  following,  definitions  and  meaning  of  terminology,  such  as 

ZV(P ),  V(P,  5),  are  as  in  the  two-dimensional  case  with  the  words  poly- 
hedron and  space  replacing  the  words  polygon  and  plane,  respectively,  when 
necessary.  However,  points  in  S2,  S3,  S4  require  some  additional  definitions 
and  terminology  which  we  present  separately  for  the  purpose  of  clarity.  Most 
importantly,  we  define  symbols  C"(P)  and  d\  for  each  point  P in  S2US3US4, 
and  describe  what  it  means  to  say  that  P ‘has  been  closed.’ 

Definitions  and  terminology  for  S2.  Let  P be  a point  in  S2.  Let  F be  the  face 
of  the  cube  closest  to  P.  Let  l be  the  ray  with  origin  P that  is  perpendicular 
to  F.  Let  P'  be  the  point  at  which  / intersects  F.  Let  m be  a line  through 
P'  that  is  perpendicular  to  an  edge  of  the  cube  in  F . Let  H"  be  the  plane 
parallel  to  F that  contains  the  point  [P  + P')/ 2.  We  define  C"{P)  as  that 
closed  half-space  determined  by  H"  that  contains  P . 

Assume  P is  the  currently  active  point,  and  P[,  i = 1,...,8,  are  points 
in  F contained  in  the  octants  around  P'  as  shown  in  Figure  6.  We  say  that 
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at  the  current  time  P has  been  closed  and  that  {P/,z  = closes  P 

if  there  exist  points  Pi,  i = 1, . . . , 8,  in  Zl  U ZV(P)  such  that  the  rays  PP;, 
i = 1, . . . , 8,  intersect  F at  the  points  P/,  i — 1, . . . , 8,  respectively.  Assuming 
P has  been  closed  we  define  d[  at  the  current  time  as  the  smallest  value  of 
dmax(P,{Pl,i  = 1,...,8})  for  {P/,z  = 1,...,8}  in  the  family  of  sets  that 
close  P. 

Definitions  and  terminology  for  S3.  Let  P be  a point  in  S 3.  Let  Pj,  j = 1,2, 
be  the  two  faces  of  the  cube  closest  to  P.  For  each  j,  j = 1,2,  let  lj  be  the 
ray  with  origin  P that  is  perpendicular  to  Fj.  For  each  j,  j = 1,2,  let  Pj  be 
the  point  at  which  intersects  Fj.  Let  Pjo,  j = 1,2,  be  the  vertices  of  the 
cube  common  to  P4  and  P2  in  the  order  shown  in  Figure  7.  Let  m be  the 
line  that  contains  the  edge  of  the  cube  common  to  Pj  and  P2.  For  each  j, 
j — 1,2,  let  m3  be  the  line  through  Pj  perpendicular  to  m.  Let  ra0  represent 
the  same  line  that  m2  represents.  For  each  j , j = 1,2,  let  Pj_i,i  and  Pj0  be 
the  closed  half-planes  determined  by  mJ_1  and  mj,  respectively,  that  contain 
Pj q.  For  each  y,  j = 1,2,  let  H'f  be  the  plane  parallel  to  Fj  that  contains 
(P  + Pj)/ 2.  We  define  C"(P)  as  the  intersection  of  the  closed  half-spaces 
determined  by  FT"  and  FFf  that  contain  P. 

Assume  P is  the  currently  active  point,  and  Pji?  j = 1,2,  z = 0, . . . , 7,  are 
points  such  that  with  Pj7  = Pj7,  for  each  j,  j = 1,2,  PjJ?  z = 1,...,6,  are 
points  in  Pj  contained  in  the  six  octants  around  Pj  as  shown  in  Figure  7,  and 
Pj  .,  7 and  Pj0  are  points  in  Pj_i,i  and  Pj0,  respectively.  We  say  that  at  the 
current  time  P has  been  closed  and  that  {Pj;,j  = 1,2 , z = 0, . . . , 7}  closes  P 
if  there  exist  points  Pj;,  j = 1,2,  i = 0, . . . , 6,  in  Zl  U ZV(P ) such  that  for 
each  j,  j = 1,2,  the  rays  PPj;,  i — 1, . . . , 6,  intersect  F3  at  the  points  Pji5  i = 

1, . . . , 6,  respectively,  and  the  ray  PPjo  intersects  Pj_i,i  and  Ejo  at  the  points 
Pj  j 7 and  Pj0,  respectively.  Assuming  P has  been  closed  we  define  dj  at  the 
current  time  as  the  smallest  value  of  dmax(P,  {P-^  j — 1,2,  i = 0,...,7}) 
for  (Pjj,  j = 1, 2,  i = 0, . . . , 7}  in  the  family  of  sets  that  close  P. 

Definitions  and  terminology  for  S4.  Let  P be  a point  in  54.  Let  Fj,  j — 1, 2,  3, 
be  the  three  faces  of  the  cube  closest  to  P in  the  order  shown  in  Figure  8. 
Let  P0  represent  the  same  face  that  P3  represents.  For  each  j,  j = 1,2,3,  let 
lj  be  the  ray  with  origin  P perpendicular  to  Pj.  For  each  j,  j = 1,2,3,  let 
Pj  be  the  point  at  which  lj  intersects  Pj.  Let  Pq  represent  the  same  point 
that  Pj  represents.  Let  Pj'  be  the  vertex  of  the  cube  common  to  P4,  P2,  and 
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Figure  6:  View  of  the  face  closest  to  a point  in  S2  that  has  been  closed. 
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Figure  7:  View  of  the  two  faces  closest  to  a point  in  53  that  has  been  closed. 
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Figure  8:  View  of  the  three  faces  closest  to  a point  in  S4  that  has  been  closed. 
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F3.  For  each  j,  j = 1,2,3,  let  rrij  be  the  line  that  contains  the  edge  of  the 
cube  common  to  Fj_l  and  Fj.  For  each  j,  j — 1,2,3,  let  rrij_1,1  and  mj0  be 
the  lines  through  Pj_1  and  Pj,  respectively,  perpendicular  to  mr  For  each  j, 
j — 1,2,3,  let  Ej_i'  1 and  Ej0  be  the  closed  half-planes  determined  by  77ij_i,i 
and  mjo,  respectively,  that  do  not  contain  P^  , and  that  are  contained  in  the 
planes  that  contain  FJ_1  and  Fj,  respectively.  For  each  j , j = 1,2,3,  let  H" 
be  the  plane  parallel  to  Fj  that  contains  (P  + Pj)/ 2.  We  define  C"{P ) as  the 
intersection  of  the  closed  half-spaces  determined  by  H",  H ",  and  H"  that 
contain  P. 

Assume  P is  the  currently  active  point,  and  Pfi,  j = 1,2,3,  i = 0,...,5, 

are  points  such  that  with  Pg5  = P^,  for  each  j , j = 1,2,3,  P-i?  i = 1,...,4, 

are  points  in  Fj  contained  in  the  four  octants  around  P-  as  shown  in  Fig- 
ure 8,  and  Pj_1 5 and  Pj0  are  points  in  Pj-1,1  and  Ej0,  respectively.  We 
say  that  at  the  current  time  P has  been  closed  and  that  {Pjiij  — 1>2,3, 
i = 0, . . . , 5}  closes  P if  there  exist  points  Pji,  j = 1, 2,  3,  i = 0, . . . , 4,  in 

U ZV{P ) such  that  for  each  j,  j = 1,2,3,  the  rays  PPji , i = 1,...,4, 

intersect  Fj  at  the  points  Pj-,  i = 1,...,4,  respectively,  and  the  ray  PPjo 
intersects  Pj_i,i  and  Pj0  at  the  points  Pj  j 5 and  P^,  respectively.  Assum- 
ing P has  been  closed  we  define  d[  at  the  current  time  as  the  smallest  value 
of  dmax(P , {P'i,j  = 1,2, 3, i = 0, ...  ,5})  for  {Pj;,  j = 1, 2, 3,  i = 0, . . . , 5}  in 
the  family  of  sets  that  close  P. 

A modified  version  of  Bowyer’s  three-dimensional  insertion  process  [2]  is 
used  in  what  follows.  It  is  the  obvious  generalization  to  three  dimensions  of 
the  modified  version  of  Bowyer’s  two-dimensional  insertion  process. 

Start  of  algorithm. 

Step  1.  Assign  points  to  cells  and  select  first  class  1 cell  to  be  activated. 

Let  M be  the  floor  of  N1^3. 

Partition  the  cube  into  M3  equal-sized  cubic  cells. 

Determine  inner  and  outer  cells. 

Assign  each  point  in  5 to  a cell. 

For  each  cell,  list  the  points  assigned  to  it. 

Determine  the  centermost  cell. 

If  the  centermost  cell  is  empty  then  go  to  Step  2. 

Else  designate  this  cell  as  the  currently  active  cell. 
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Go  to  Step  3. 

Step  2.  Select  next  inner  ( class  1,  2,  3,  or  4)  cell. 

If  all  inner  cells  have  been  activated  then  go  to  Step  8. 

Else  choose  the  next  inner  cell  to  be  activated. 

If  this  cell  is  empty  then  go  to  Step  2. 

Else  designate  this  cell  as  the  currently  active  cell. 

Determine  class  for  currently  active  cell. 

If  class  1 then  go  to  Step  3. 

Else  go  to  Step  4. 

Step  3.  Construct  Voronoi  polyhedron  of  a point  in  S\. 

Let  P be  a point  assigned  to  the  currently  active  cell  that  has  not  been 
processed. 

Designate  P as  the  currently  active  point. 

Start  spiral  search  around  P and  construct  V1. 

Update  V1  and  each  subsequent  V1  as  appropriate. 

For  each  \n  compute  Dt  = dmax(P.Vt)  and  dt. 

Terminate  search  when  one  of  the  following  criteria  is  met. 

1.  2 Dt  < dt. 

2.  All  cells  in  the  cube  have  been  searched. 

Upon  termination  go  to  Step  7. 

Step  4.  Begin  construction  of  Voronoi  polyhedron  of  a point  in  S n,  S 3,  or 
£4. 

Let  P be  a point  assigned  to  the  currently  active  cell  that  has  not  been 
processed. 

Designate  P as  the  currently  active  point. 

Start  spiral  search  around  P and  construct  V1 . 

Update  V1  and  each  subsequent  Vl  as  appropriate. 

For  each  \rt  compute  Dt  = dmax(P,Vt ) and  dt. 

Terminate  search  when  one  of  the  following  criteria  is  met. 

1.  2 Dt  < dt. 

2.  All  cells  in  the  cube  have  been  searched. 
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3.  P has  been  closed. 


Upon  termination,  if  neither  criterion  1 nor  criterion  2 has  been  met 
then  go  to  Step  5. 

Else  go  to  Step  7. 

Step  5.  Continue  construction  of  Voronoi  polyhedron  of  a point  in  S 2,  S3, 
or  S4. 

Resume  spiral  search  around  P. 

Update  each  V1  as  appropriate. 

For  each  V1  compute  Dt  = dmax(P,\n),  dt,  and  d[. 

Terminate  search  when  one  of  the  following  criteria  is  met. 

1.  2 Dt  < dt . 

2.  All  cells  in  the  cube  have  been  searched. 

3.  y/2  d[  < dt. 

Upon  termination,  if  neither  criterion  1 nor  criterion  2 has  been  met 
then  go  to  Step  6. 

Else  go  to  Step  7. 

Step  6.  Complete  construction  of  Voronoi  polyhedron  of  a point  in  So,  S3, 
or  S 4. 

Determine  C"(P). 

Resume  spiral  search  around  P. 

Update  each  V1  as  appropriate. 

For  each  \n  compute  Dt  = dmax(P.Vt  D C"(P))  and  dt. 

Terminate  search  when  one  of  the  following  criteria  is  met. 

1.  2 Dt  < dt. 

2.  All  cells  in  the  cube  have  been  searched. 

Upon  termination  go  to  Step  7. 

Step  7.  Save  Voronoi  polyhedron  of  a point  assigned  to  an  inner  cell. 
Identify  U(F,  S)  with  Vt. 

Mark  P as  processed  and  save  V(P,  S). 

For  each  P'  in  U(P,  S)  that  has  not  been  processed  let 
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ZV{P')  = Zv{P')u{P}. 

Determine  whether  currently  active  cell  has  been  activated. 

If  activated  then  go  to  Step  2. 

Else  if  P is  in  Si  then  go  to  Step  3. 

Else  go  to  Step  4. 

Step  8.  Construct  and  save  Voronoi  polyhedra  of  points  in  S 5. 

Determine  S5. 

If  S5  is  empty  then  stop. 

Else  let  Z5  = U PesBZv(P). 

Perform  insertion  process  on  S5  U Z5. 

For  each  P in  S5  identify  V(P,S)  with  Vr(P,  S5  U Z5). 

For  each  P in  S5  mark  P as  processed  and  save  V(P,  5). 

Stop. 

End  of  algorithm. 

Justification  of  algorithm.  Because  of  similarities  with  the  two  dimensional 
method,  we  need  not  justify  that  the  above  algorithm  constructs  the  Voronoi 
polyhedra  of  points  in  Si  or  S5.  In  addition,  because  of  similarities  among 
52,  S3,  and  54,  we  only  justify  that  it  constructs  the  Voronoi  polyhedra  of 
points  in  5 2. 

Let  F , l,  P' , m,  H" , C”{P ) be  as  defined  above  for  a point  P in  5 2.  Assume 
that  P is  the  currently  active  point  and  that  it  has  been  closed.  In  addition, 
assume  that  at  the  current  time  y/2  d[  < dt  and  that  a set  {P/,z  = 1, . . . , 8} 
of  points  in  F closes  P with  d[  = dmax(P , {P/,  i = 1, . . . , 8}).  Let  V be  that 
part  of  V1  that  is  not  contained  in  C"(P).  Assume  S \ Zl  U ZV[P ) U {P}  is 
not  empty  and  O'  is  a point  in  this  set.  Let  H'  be  the  plane  that  perpen- 
dicularly bisects  O' P , and  let  C'  be  the  open  half-space  determined  by  H' 
that  contains  P.  We  show  that  C'  contains  V,  so  that  O'  does  not  affect  V , 
and  thus  the  termination  criteria  may  change  from  those  in  Step  5 to  those 
in  Step  6. 

Let  Pi,  i = 1 , . . . , 8,  be  points  in  Z 1 U ZV(P)  such  that  the  rays  PP,, 
i = 1,...,8,  intersect  P at  the  points  P-,  i = 1,...,8,  respectively.  We 
assume  P/  = Pt  for  each  i,i  — 1, . . . , 8,  since,  at  worst,  V is  just  a larger  set. 
We  assume  P/  ^ P' , i = 1,. . . ,8,  since,  otherwise,  V equals  the  empty  set. 
Finally,  we  assume  P is  not  a convex  combination  of  P'  and  O'. 
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We  first  argue  that  we  may  assume  that  O'  is  contained  in  F . For  suppose 
O'  is  not  in  F . Let  Q'0  be  the  point  in  the  plane  that  contains  F such  that 
O'  and  Q'0  are  equidistant  from  P,  and  Q'q  lies  in  the  half-plane  that  contains 
Q'  and  which  is  determined  by  a line  through  P perpendicular  to  F.  Such  a 
O'0  exists  because  the  perpendicular  distance  from  P to  F is  less  than  d[  and 
dist(Q' , P)  > y/2  d't.  Let  H'Q  be  the  plane  that  is  the  perpendicular  bisector 
of  Q'0P.  Let  C'q  be  the  open  half-space  determined  by  H'Q  that  contains  P . 
It  follows,  as  shown  in  Figure  9,  that  if  C'Q  contains  V,  which  lies  above  H'\ 
then  so  does  C'.  Thus,  without  any  loss  of  generality,  we  assume  that  Q' 
is  contained  in  F.  Further,  we  assume  that  P' P[  / P' P2  and  that  O'  is 
contained  in  the  region  in  F obtained  by  a clockwise  rotation  along  F from 
P'  P2  to  P' P[  (refer  back  to  Figure  6). 

Let  H[  and  H2  be  the  planes  that  are  the  perpendicular  bisectors  of  P[P 
and  P2P , respectively.  Let  R be  the  region  which  is  the  intersection  of  the 
closed  half-space  determined  by  H"  that  does  not  contain  P and  the  closed 
half-spaces  determined  by  H[  and  H2  that  contain  P.  Clearly  R contains  V. 
We  will  show  that  C'  contains  V by  showing  that  R is  the  convex  hull  of  a 
region  (A"")  and  a ray  ( r "),  both  of  which  lie  in  C' . Since  C’  is  convex,  C' 
then  contains  R and  the  result  follows. 

The  following  definitions  will  be  useful.  Let  H"’  be  the  plane  that  contains  P 
and  is  parallel  to  H" . Let  Q"\  P"\  and  P2"  be  the  perpendicular  projections 
onto  H"’  of  O',  Pj,  and  P2',  respectively.  Let  d'0 , d'1?  d2,  d'^' , d"\  and  d'"  be  the 
values  of  dist(Q' , P),  dist(P(,  P),  dist(P2,  P),  dist(Q'" , P),  dist(P"' , P),  and 
dist(P"' , P),  respectively.  Let  /q,  and  l"  be  the  lines  that  are  the  inter- 
sections of  H"  with  if',  H[,  and  id2,  respectively.  Let  and  V2  be  the 

lines  in  H'"  that  perpendicularly  bisect  0'"P,  P{" P , and  P"'P<,  respectively. 
Let  P"  be  the  perpendicular  projection  of  P onto  H" . The  region  K"  is 
defined  to  be  the  wedge  obtained  by  intersecting  the  half-planes  in  H"  de- 
termined by  l"  and  l"  that  contain  P"  (Figure  10).  Another  region,  K"\  is 
similarly  defined  with  respect  to  iP",  /2",  and  P (Figure  11). 

We  begin  the  task  of  showing  that  C'  contains  the  region  K"  by  arguing 
that  dq  > y/2  dmax(P,  {P"' ■>  P”'})- 

(O'  + P)/2,  (P[  + P)/2,  (P'  + P)/2,  (O'"  + P)/2,  (PI"  + P)/2,  (P'"  + P)/2 
are  contained  in  1'0',  l ",  /2,  /q",  V" , respectively,  so  that  from  the  definitions 

of  0'",  P"',  and  P2",  it  follows  by  similar  triangles  that  and  l2  are  the 

perpendicular  projections  onto  PL'"  of  /J,  /",  and  /",  respectively.  If  h'  is  the 
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Figure  9:  Perpendicular  view  of  unique  plane  that  contains  P,  Q'0,  and  Q' . 
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/ 

/ 
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Figure  10:  Perpendicular  view  of  H" . The  shaded  area  is  the  wedge  K . 
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p'"  / 

/ 

/ 


Figure  11:  Perpendicular  view  of  H"'.  The  shaded  area  is  the  wedge  K"' . 
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value  of  dist{P\  P)  then,  since  d'0  > \/2d't,  we  have 

«T)2  + (*')2  = K)2 
> 2 «)2 

> 2 K)2 

= 2(Ki")2  + W2) 

> 2 (d'?y  + (h')2 . 

Thus,  d"'  > y/2  d'" . Similarly  d g"  > \/2d2",  which  shows 

d"'  > \/2  dmax(P:  {P”' ■>  P"'})- 

Since  Pf  and  P2  are  in  contiguous  octants,  it  follows  that  I'f  does  not  intersect 
K1".  In  addition,  since  K’"  is  the  perpendicular  projection  onto  H1"  of  if", 
it  also  follows  that  /J  does  not  intersect  K"  (Figures  10  and  11).  In  other 
words,  by  the  definition  of  C'  contains  K". 

Next,  we  define  r"  and  show  it  is  contained  in  C’ . 

Let  H*  be  the  unique  plane  that  contains  P,  Pj,  and  P2.  Let  C*  be  the 
closed  half-space  determined  by  H*  that  contains  PL  Since,  in  particular, 
dj,  > dmax(P,  {P[,  P2})?  if  follows  that  some  convex  combination,  say  T',  of 
P[  and  P2  is  also  a convex  combination  of  P'  and  O'  with  T'  ^ O'  (Figure  12). 
Thus,  since  P'  is  an  interior  point  of  C*,  and  T'  is  contained  in  if*,  it  follows 
that  O'  does  not  belong  to  C*. 

Let  / " be  the  line  that  is  the  intersection  of  the  planes  H[  and  if2.  Let  H** 
be  the  plane  that  contains  P and  O'  and  that  is  perpendicular  to  H* . Let 
/**  be  the  perpendicular  projection  onto  H**  of  l".  Since  l"  is  perpendicular 
to  H*  (Figure  13),  so  is  /**.  Thus,  as  shown  in  Figure  14,  /**  contains  a ray 
r**  that  lies  wholly  in  C'  D C* . Therefore,  by  the  definition  of  /’*,  it  follows 
that  l"  contains  a ray  r"  that  is  also  contained  in  C'  D C* . 

Since  both  K"  and  r"  lie  in  C',  their  convex  hull  R lies  in  the  convex  region 
C',  and  our  assertion  is  proved. 

Proof  of  complexity.  We  now  assume  that  the  points  in  5 have  been  cho- 
sen from  a uniform  distribution  in  the  cube.  Under  this  assumption,  we 
claim  that  the  expected  time  involved  with  the  above  algorithm  is  O(N) 
and  0(N 4^3)  for  obtaining  all  Voronoi  polyhedra  of  points  in  5 \ 5s  and  S5, 
respectively.  Again,  because  of  similarities  with  the  two-dimensional  algo- 
rithm, we  only  present  the  proof  for  S5. 
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Figure  12:  Perpendicular  view  of  face  of  cube  closest  to  P.  P{P^  and  P'Q' 
intersect  at  the  point  T1. 
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p; 


Figure  13:  Perpendicular  view  of  plane  H* . /"  is  perpendicular  to  H * through 
the  point  T* . 
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Figure  14:  Perpendicular  view  of  plane  H** . The  shaded  area  is  C'  H C*. 
The  ray  r*m  is  the  solid  portion  of  the  line  /**. 
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Let  Z5  be  as  defined  in  Step  8 of  the  algorithm.  We  first  show  that  the 
expected  number  of  points  in  55  U Z5  is  0(N2! 3). 

Clearly,  the  expected  number  of  points  in  55  is  0(X2/3).  Thus,  it  suffices  to 
prove  that  the  expected  number  of  points  in  Z5  is  0(N 2//3). 

In  what  follows,  given  G,  a finite  nonempty  set  in  Euclidean  space,  and 
Gi,  G2,  nonempty  subsets  of  G,  we  define  X(G,  Gi,G2)  as  the  number  of 
points  in  G2  that  are  Voronoi  neighbors  in  G of  points  in  Gi.  Accordingly, 
E(N(G , Gi,  G2))  is  defined  as  the  expected  value  of  X(G,  Gi,  G2). 

The  following  observation  is  crucial  to  our  proof: 

Let  X,  Y be  finite  nonempty  sets  in  Euclidean  space.  Let  X'  be  a nonempty 
subset  of  X . Then 

X(X,X',X)  < N(X  U V,X',X)  + N(X  U Y,Y,X). 

We  prove  a stronger  result.  Let  P and  P'  be  points  in  X and  X',  respectively. 
We  show  that  if  P'  is  a Voronoi  neighbor  in  Ar  of  P then  either  P'  is  a Voronoi 
neighbor  in  A"  U Y of  P or  there  exists  at  least  one  point  in  Y which  is  a 
Voronoi  neighbor  in  A"  U Y of  P.  If  not,  all  Voronoi  neighbors  in  X U Y of  P 
lie  in  X,  none  of  which  is  P\  But  this  implies  that  the  Voronoi  polyhedron 
of  P relative  to  X coincides  with  the  Voronoi  polyhedron  of  P relative  to 
X U V,  so  that  P'  cannot  be  a Voronoi  neighbor  in  A"  of  P,  a contradiction. 
Let  B be  the  cube  from  which  the  set  S has  been  chosen.  Let  B'  be  the  cube 
obtained  by  surrounding  B with  L(N ) + 1 additional  layers  of  cells.  Define 
cells  in  B'  \ B as  exterior  cells.  Let  Z be  the  set  whose  members  are  the 
centroids  of  the  exterior  cells. 

It  follows  from  the  observation  above  that 


N{S,  5S,  S)  < N(S  U Z,  S5,  S)  + N(S  U Z,  Z,  5). 

Thus,  since  each  term  above  is  a random  variable  that  depends  solely  on  the 
choice  of  5,  we  must  have 

E(N(S,  S5, 5))  < E{N(S  U Z,  S5,  S))  + E{N(S  U Z,  Z,  5)). 

Let  Z'  be  the  set  of  points  in  Z contained  in  the  first  layer  of  exterior  cells 
that  surrounds  B.  From  the  definition  of  Z,  points  in  Z'  are  the  only  points 
in  Z that  can  be  Voronoi  neighbors  in  5 U Z of  points  in  5.  Clearly,  points 
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in  S5  U Z ' are  not  contained  in  any  of  the  outermost  L(N ) layers  of  cells  of 
B\  so  that  the  expected  number  of  Voronoi  neighbors  in  5 U Z of  a point  in 
S5  U Z ' is  0(1)  [l].  Since  the  expected  number  of  points  in  55  U Z'  is  0(iV2//3), 
this  implies 

E(N{S  u Z,  S5, 5))  = E(N(S  U Z,Z,S))  = 0{ N2'3), 

so  that 

E(N(S,  Ss,  5))  < 0(N2/1)  + 0(N2/3)  = 0(N2/3). 

Since  E(N(S , S5,  5))  is  by  definition  an  upper  bound  on  the  expected  number 
of  points  in  Z 5,  our  assertion  follows. 

Finally,  we  show  that  the  expected  time  required  to  perform  the  insertion 
process  of  Step  8 is  bounded  above  by  0(iV4//3). 

Assume,  inductively,  that  the  Voronoi  diagram  for  k — 1 points  in  S 5 U has 

been  constructed.  Since  the  expected  number  of  points  in  S5UZ5  is  0(V2/3), 
it  follows  that  0(N 2/3)  is  an  upper  bound  on  the  expected  time  required  to 
find  that  one  of  the  k — 1 points  closest  to  an  additional  kth  point,  and  that 
the  expected  number  of  vertices  of  any  polyhedron  in  the  Voronoi  diagram  for 
the  k — 1 points  is  at  most  0(  V2y/3).  Thus,  the  expected  time  required  to  find 
a vertex  of  the  Voronoi  diagram  for  the  k — 1 points  that  does  not  belong 
to  the  Voronoi  diagram  for  all  k points  is  at  most  0(vV2/3).  This  implies 
that  the  expected  time  required  by  the  first  step  of  the  insertion  process,  as 
described  in  the  introduction,  is  0(N 2//3)  for  each  point  in  S5  U Z5,  and 

0(iV2/3).0(iV2/3)  = 0(iV4/3) 

for  all  of  55  U Z5. 

Next,  we  show  that  the  expected  time  required  by  the  remainder  of  the 
insertion  process  is  also  0(N 4^3)  for  all  of  S5  U Z5.  As  described  in  the 
introduction,  and  again  assuming  that  the  Voronoi  diagram  for  k — 1 points 
in  S5UZ5  has  been  constructed,  the  remainder  of  the  insertion  process  consists 
of  deleting  those  vertices  in  the  Voronoi  diagram  for  the  k — 1 points  that  lie 
in  the  Voronoi  polyhedron  of  the  kth  point,  and  adding  to  the  diagram  the 
vertices  of  the  Voronoi  polyhedron  of  the  kth  point.  Let  r be  the  number  of 
points  in  S5  U and  for  each  fc,  1 < k < r,  let  d(k ) and  a(k)  be  the  number 
of  vertices  that  are  deleted  and  added,  respectively,  when  updating  for  the 
kth  point  in  S5  U Z5  the  Voronoi  diagram  for  the  first  k — 1 points.  Clearly, 
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for  each  k,  1 < k < r,  the  expected  value  of  a(k)  is  at  most  0(iV2^3),  so 
that  the  expected  value  of  Y7k=i  aW  mosf  0(N 4)/3).  Also,  since  deleted 
vertices  come  from  the  set  of  previously  added  vertices,  we  must  have 

Y.d(k)  < Y,a(k), 

k— 1 k — 1 

and  the  desired  result  follows. 


5.  Computational  experience 

The  two-dimensional  algorithm  was  implemented  in  standard  FORTRAN 
on  a Control  Data  Cyber  2051  at  the  National  Bureau  of  Standards.  The 
storage  of  the  data  structure  and  the  treatment  of  degenerate  vertices  were 
based  on  considerations  from  [2].  Table  1 shows  the  computing  times  per 


N 

Proposed  algorithm  Bowyer’s  algorithm 

8,100 

1.27 

X 

1(T3 

1.07 

X 

10"3 

14,400 

1.37 

X 

10~3 

1.29 

X 

10'3 

52,900 

1.30 

X 

10~3 

1.97 

X 

10"3 

102,400 

1.28 

X 

10"3 

2.59 

X 

10~3 

122,500 

1.28 

X 

10~3 

2.80 

X 

10"3 

Table  1:  Computing  time  per  Voronoi  polygon. 


Voronoi  polygon  for  the  algorithm  when  applied  to  several  regular  nontrivial 
configurations.  The  times  obtained  when  applying  the  modified  version  of 
Bowyer’s  algorithm  to  the  same  sets  are  also  shown.  The  unit  of  computing 
time  is  given  in  CPU  seconds. 
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